//#define ALLOW_FLOAT128 1
#ifndef __STDC_WANT_IEC_60559_TYPES_EXT__
#define __STDC_WANT_IEC_60559_TYPES_EXT__ 1
#endif

#include <cfloat>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <ctgmath>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <limits>
#ifdef __has_include
# if __has_include(<stdfloat>)
//#   include <stdfloat>
# endif
#endif
#if defined(FLT128_MIN) && USE_QUADMATH
#include <quadmath.h>
#endif
#include "arm/neon.h"

using namespace std;

#define FLAGS_NONE 0
#define FLAGS_LAST 1

#define BUFSIZE (1024 * 1024)

typedef void (*BenchmarkFunction)(const void* pbuf, size_t len, int flags);

static void* pbuf = NULL;
static volatile float dstFloat = 0;
static volatile double dstDouble = 0;
static volatile long double dstLongDouble = 0;
static volatile double mopsBaseline = 0;

static const char hex_table[] = "0123456789ABCDEF";

template<typename T>
size_t toHex(char* pbuf, const T* pvalue) {
    const size_t count = sizeof(T);
    const size_t len = sizeof(T) * 2;
    //memset(pbuf, '0', len);
    const char* p = (const char *)pvalue;
#if Q_BYTE_ORDER == Q_LITTLE_ENDIAN
    char* q = pbuf + (count-1)*2;
    int cchItem = -2;
#else
    char* q = pbuf;
    int cchItem = 2;
#endif
    for(size_t i =0; i<count; ++i) {
        int n = ((int)*p) & 0xFF;
        q[0] = hex_table[n >> 4];
        q[1] = hex_table[n & 0xF];
        // Next.
        ++p;
        q += cchItem;
    }
    return len;
}

template<typename T>
ostream& outOnlyHex(ostream& os, const T& value) {
    size_t len = sizeof(T) * 2;
    char buf[len+1];
    //memset(buf, 0, sizeof(buf));
    size_t n = toHex(buf, &value);
    buf[n] = '\0';
    os << "0x";
    os << buf;
    //os << "";
    return os;
}

template<typename T>
ostream& outWithHex(ostream& os, T value) {
    os << value << "(";
    outOnlyHex(os, value);
    os << ")";
    return os;
}

#ifdef QUADMATH_H
ostream& outValue(ostream& os, __float128 value) {
    const size_t size = 50;
    char buf[size+1];
    int n = quadmath_snprintf(buf, size, "%.33Qg", value);
    buf[n] = '\0';
    os << buf;
    return os;
}

ostream& outWithHex(ostream& os, __float128 value) {
    outValue(os, value) << "(";
    outOnlyHex(os, value);
    os << ")";
    return os;
}
#endif // QUADMATH_H

double benchmarkCall(double* ms, const char* name, BenchmarkFunction func, void* pbuf, size_t len) {
    const int count = 1000;
    clock_t clockBegin, clockDiff;
    func(pbuf, len, FLAGS_NONE);
    // body.
    clockBegin = clock();
    for(int i=0; i<count; ++i) {
        func(pbuf, len, FLAGS_NONE);
    }
    clockDiff = clock() - clockBegin;
    *ms = (double)((long double)clockDiff * 1000 / (long double)(CLOCKS_PER_SEC));
    double mops = (double)( (long double)(len * count)/(1024*1024) / ((long double)clockDiff / CLOCKS_PER_SEC) );
    // done.
    func(pbuf, len, FLAGS_LAST);
    return mops;
}

void benchmarkCallShow(const char* name, BenchmarkFunction func, void* pbuf, size_t len) {
    double mops, ms;
    mops = benchmarkCall(&ms, name, func, pbuf, len);
    //cout << fixed;
    cout << name << ": " << ms <<  "ms, " << mops << " mops";
    //cout << scientific;
    if (mopsBaseline > 0) {
        double scale = mops / mopsBaseline;
        cout << ", scale=" << scale;
    } else {
        mopsBaseline = mops;
    }
}

#ifdef __STDCPP_FLOAT16_T__
static volatile float16_t dstFloat16 = 0;

void sumFloat16(const void* pbuf, size_t len, int flags) {
    float16_t rt = 0;
    const float16_t* p = (const float16_t*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstFloat16 = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumFloat16OMP(const void* pbuf, size_t len, int flags) {
    float16_t rt = 0;
    const float16_t* p = (const float16_t*)pbuf;
    size_t i;
#pragma omp parallel for reduction(+ : rt)
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstFloat16 = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumFloat16Simd(const void* pbuf, size_t len, int flags) {
    float16_t rt = 0;
    const size_t vectorWidth = sizeof(simde_float16x8_t)/sizeof(float16_t);
    const float16_t* p = (const float16_t*)pbuf;
    const float16_t* pend = p + len;
    const float16_t* pendv = p + (len/vectorWidth)*vectorWidth;
    size_t i;
    simde_float16x8_t vrt = simde_vdupq_n_f16(0);
    while(p < pendv) {
        //rt += *p * *p;
        simde_float16x8_t va = simde_vld1q_f16(p);
        simde_float16x8_t vc = simde_vmulq_f16(va, va);
        vrt = simde_vaddq_f16(vrt, vc);
        p += vectorWidth;
    }
    while(p < pend) {
        rt += *p * *p;
        ++p;
    }
    //rt += sum(vrt);
    p = (const float16_t*)(const void*)&vrt;
    for(i=0; i < vectorWidth; ++i) {
        rt += *p++;
    }
    // done.
    dstFloat16 = rt;
    if (flags & FLAGS_LAST) {
    }

}

void sumFloat16To32(const void* pbuf, size_t len, int flags) {
    float rt = 0;
    const float16_t* p = (const float16_t*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt += (float)(*p) * *p;
        ++p;
    }
    // done.
    dstFloat = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumFloat16To32OMP(const void* pbuf, size_t len, int flags) {
    float rt = 0;
    const float16_t* p = (const float16_t*)pbuf;
    size_t i;
#pragma omp parallel for reduction(+ : rt)
    for(i=0; i < len; ++i) {
        rt += (float)(*p) * *p;
        ++p;
    }
    // done.
    dstFloat = rt;
    if (flags & FLAGS_LAST) {
    }
}

float addCross_f32(simde_float32x4_t a) {
    float rt = 0;
    const size_t vectorWidth = 4;//sizeof(simde_float32x4_t)/sizeof(float);
    const float* p = (const float*)(const void*)&a;
    size_t i;
    for(i=0; i < vectorWidth; ++i) {
        rt += *p++;
    }
    return rt;
}

void sumFloat16To32Simd(const void* pbuf, size_t len, int flags) {
    float rt = 0;
    const size_t vectorWidth = 8;//sizeof(simde_float16x8_t)/sizeof(float16_t);
    const float16_t* p = (const float16_t*)pbuf;
    const float16_t* pend = p + len;
    const float16_t* pendv = p + (len/vectorWidth)*vectorWidth;
    size_t i;
    simde_float32x4_t vrt = simde_vdupq_n_f32(0.0f);
    simde_float32x4_t vrt1 = simde_vdupq_n_f32(0.0f);
    while(p < pendv) {
        //rt += *p * *p;
        simde_float16x8_t va = simde_vld1q_f16(p);
        simde_float32x4_t va0 = simde_vcvt_f32_f16(simde_vget_low_f16(va));
        simde_float32x4_t va1 = simde_vcvt_high_f32_f16(va);
        simde_float32x4_t vc0 = simde_vmulq_f32(va0, va0);
        simde_float32x4_t vc1 = simde_vmulq_f32(va1, va1);
        vrt = simde_vaddq_f32(vrt, vc0);
        vrt1 = simde_vaddq_f32(vrt1, vc1);
        p += vectorWidth;
    }
    vrt = simde_vaddq_f32(vrt, vrt1);
    while(p < pend) {
        rt += (float)(*p) * *p;
        ++p;
    }
    rt = addCross_f32(vrt);
    // done.
    dstFloat = rt;
    if (flags & FLAGS_LAST) {
    }

}

void callFloat16(void* pbuf, size_t len) {
    float16_t* p = (float16_t*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        *p++ = (float16_t)((rand() & 0xFFFFU) / 65536.0f);
    }
    mopsBaseline = 0;
    // run.
    benchmarkCallShow("sumFloat16", sumFloat16, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat16) << endl;
    benchmarkCallShow("sumFloat16OMP", sumFloat16OMP, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat16) << endl;
    benchmarkCallShow("sumFloat16Simd", sumFloat16Simd, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat16) << endl;
    benchmarkCallShow("sumFloat16To32", sumFloat16To32, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat) << endl;
    benchmarkCallShow("sumFloat16To32OMP", sumFloat16To32OMP, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat) << endl;
    benchmarkCallShow("sumFloat16To32Simd", sumFloat16To32Simd, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat) << endl;
}

#endif // __STDCPP_FLOAT16_T__

#ifdef __STDCPP_FLOAT16_T__
static volatile bfloat16_t dstBFloat16 = 0;

void sumBFloat16(const void* pbuf, size_t len, int flags) {
    bfloat16_t rt = 0;
    const bfloat16_t* p = (const bfloat16_t*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstBFloat16 = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumBFloat16OMP(const void* pbuf, size_t len, int flags) {
    bfloat16_t rt = 0;
    const bfloat16_t* p = (const bfloat16_t*)pbuf;
    size_t i;
#pragma omp parallel for reduction(+ : rt)
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstBFloat16 = rt;
    if (flags & FLAGS_LAST) {
    }
}

void callBFloat16(void* pbuf, size_t len) {
    bfloat16_t* p = (bfloat16_t*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        *p++ = (bfloat16_t)((rand() & 0xFFFFU) / 65536.0f);
    }
    mopsBaseline = 0;
    // run.
    benchmarkCallShow("sumBFloat16", sumBFloat16, pbuf, len); cout << " // result="; outWithHex(cout, dstBFloat16) << endl;
    benchmarkCallShow("sumBFloat16OMP", sumBFloat16OMP, pbuf, len); cout << " // result="; outWithHex(cout, dstBFloat16) << endl;
}

#endif // __STDCPP_FLOAT16_T__

void sumFloat(const void* pbuf, size_t len, int flags) {
    float rt = 0;
    const float* p = (const float*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstFloat = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumFloatFma(const void* pbuf, size_t len, int flags) {
    float rt = 0;
    const float* p = (const float*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt = fmaf(*p, *p, rt);
        ++p;
    }
    // done.
    dstFloat = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumFloatOMP(const void* pbuf, size_t len, int flags) {
    float rt = 0;
    const float* p = (const float*)pbuf;
    size_t i;
#pragma omp parallel for reduction(+ : rt)
    for(i=0; i < len; ++i) {
        rt = fma(*p, *p, rt);
        ++p;
    }
    // done.
    dstFloat = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumFloatSimd(const void* pbuf, size_t len, int flags) {
    float rt = 0;
    const size_t vectorWidth = sizeof(simde_float32x4_t)/sizeof(float);
    const float* p = (const float*)pbuf;
    const float* pend = p + len;
    const float* pendv = p + (len/vectorWidth)*vectorWidth;
    size_t i;
    simde_float32x4_t vrt = simde_vdupq_n_f32(0.0f);
    while(p < pendv) {
        //rt += *p * *p;
        simde_float32x4_t va = simde_vld1q_f32(p);
        simde_float32x4_t vc = simde_vmulq_f32(va, va);
        vrt = simde_vaddq_f32(vrt, vc);
        p += vectorWidth;
    }
    while(p < pend) {
        rt += *p * *p;
        ++p;
    }
    //rt += sum(vrt);
    p = (const float*)(const void*)&vrt;
    for(i=0; i < vectorWidth; ++i) {
        rt += *p++;
    }
    // done.
    dstFloat = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumDouble(const void* pbuf, size_t len, int flags) {
    double rt = 0;
    const double* p = (const double*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstDouble = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumDoubleFma(const void* pbuf, size_t len, int flags) {
    double rt = 0;
    const double* p = (const double*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt = fma(*p, *p, rt);
        ++p;
    }
    // done.
    dstDouble = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumDoubleOMP(const void* pbuf, size_t len, int flags) {
    double rt = 0;
    const double* p = (const double*)pbuf;
    size_t i;
#pragma omp parallel for reduction(+ : rt)
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstDouble = rt;
    if (flags & FLAGS_LAST) {
    }
}
void sumLongDouble(const void* pbuf, size_t len, int flags) {
    long double rt = 0;
    const long double* p = (const long double*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstLongDouble = rt;
    if (flags & FLAGS_LAST) {
        //cout << "sumLongDouble result: " << dstLongDouble << endl;
    }
}

void sumLongDoubleFma(const void* pbuf, size_t len, int flags) {
    long double rt = 0;
    const long double* p = (const long double*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt = fmal(*p, *p, rt);
        ++p;
    }
    // done.
    dstLongDouble = rt;
    if (flags & FLAGS_LAST) {
        //cout << "sumLongDouble result: " << dstLongDouble << endl;
    }
}

void sumLongDoubleOMP(const void* pbuf, size_t len, int flags) {
    long double rt = 0;
    const long double* p = (const long double*)pbuf;
    size_t i;
#pragma omp parallel for reduction(+ : rt)
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstLongDouble = rt;
    if (flags & FLAGS_LAST) {
        //cout << "sumLongDouble result: " << dstLongDouble << endl;
    }
}

void callFloat(void* pbuf, size_t len) {
    float* p = (float*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        *p++ = (float)(rand() & 0xFFFFFFU) / 65536.0f;
    }
    mopsBaseline = 0;
    // run.
    benchmarkCallShow("sumFloat", sumFloat, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat) << endl;
    benchmarkCallShow("sumFloatFma", sumFloatFma, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat) << endl;
    benchmarkCallShow("sumFloatOMP", sumFloatOMP, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat) << endl;
    benchmarkCallShow("sumFloatSimd", sumFloatSimd, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat) << endl;
}

void callDouble(void* pbuf, size_t len) {
    double* p = (double*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        *p++ = (double)(rand() & 0xFFFFFFU) / 65536.0;
    }
    mopsBaseline = 0;
    // run.
    benchmarkCallShow("sumDouble", sumDouble, pbuf, len); cout << " // result="; outWithHex(cout, dstDouble) << endl;
    benchmarkCallShow("sumDoubleFma", sumDoubleFma, pbuf, len); cout << " // result="; outWithHex(cout, dstDouble) << endl;
    benchmarkCallShow("sumDoubleOMP", sumDoubleOMP, pbuf, len); cout << " // result="; outWithHex(cout, dstDouble) << endl;
}

void callLongDouble(void* pbuf, size_t len) {
    long double* p = (long double*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        *p++ = (long double)(rand() & 0xFFFFFFU) / 65536.0;
    }
    mopsBaseline = 0;
    // run.
    benchmarkCallShow("sumLongDouble", sumLongDouble, pbuf, len); cout << " // result="; outWithHex(cout, dstLongDouble) << endl;
    //sumLongDoubleFma: 19144ms, 52.2357 mops
    //benchmarkCallShow("sumLongDoubleFma", sumLongDoubleFma, pbuf, len); cout << " // result="; outWithHex(cout, dstLongDouble) << endl;
    benchmarkCallShow("sumLongDoubleOMP", sumLongDoubleOMP, pbuf, len); cout << " // result="; outWithHex(cout, dstLongDouble) << endl;
}

#ifdef QUADMATH_H
static volatile __float128 dstFloat128 = 0;

void sumFloat128(const void* pbuf, size_t len, int flags) {
    __float128 rt = 0;
    const __float128* p = (const __float128*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstFloat128 = rt;
    if (flags & FLAGS_LAST) {
    }
}

void sumFloat128OMP(const void* pbuf, size_t len, int flags) {
    __float128 rt = 0;
    const __float128* p = (const __float128*)pbuf;
    size_t i;
#pragma omp parallel for reduction(+ : rt)
    for(i=0; i < len; ++i) {
        rt += *p * *p;
        ++p;
    }
    // done.
    dstFloat128 = rt;
    if (flags & FLAGS_LAST) {
    }
}

void callFloat128(void* pbuf, size_t len) {
    __float128* p = (__float128*)pbuf;
    size_t i;
    for(i=0; i < len; ++i) {
        *p++ = (__float128)(rand() & 0xFFFFFFU) / 65536.0f;
    }
    mopsBaseline = 0;
    // run.
    benchmarkCallShow("sumFloat128", sumFloat128, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat128) << endl;
    benchmarkCallShow("sumFloat128OMP", sumFloat128OMP, pbuf, len); cout << " // result="; outWithHex(cout, dstFloat128) << endl;
}

#endif // QUADMATH_H

void callRun() {
    pbuf = malloc(BUFSIZE * 32);



#ifdef __STDCPP_FLOAT16_T__
    callFloat16(pbuf, BUFSIZE);
#endif // __STDCPP_FLOAT16_T__
#ifdef __STDCPP_BFLOAT16_T__
    callBFloat16(pbuf, BUFSIZE);
#endif // __STDCPP_FLOAT16_T__
    callFloat(pbuf, BUFSIZE);
    callDouble(pbuf, BUFSIZE);
    callLongDouble(pbuf, BUFSIZE);
#ifdef QUADMATH_H
    callFloat128(pbuf, BUFSIZE);
#endif // QUADMATH_H

    free(pbuf);
}

int main(int argc, char *argv[])
{
    cout << "benchmarkfloat" << endl;
    cout << endl;
    srand(1);

    cout << "DECIMAL_DIG: " << DECIMAL_DIG << endl;
    cout << "FLT_RADIX: " << FLT_RADIX << endl;
    cout << "FLT_ROUNDS: " << FLT_ROUNDS << endl;
    cout << "FLT_EVAL_METHOD: " << FLT_EVAL_METHOD << endl;
#ifdef FP_FAST_FMAF
    cout << "FP_FAST_FMAF: " << FP_FAST_FMAF << endl;
#endif
#ifdef FP_FAST_FMA
    cout << "FP_FAST_FMA: " << FP_FAST_FMA << endl;
#endif
#ifdef FP_FAST_FMAL
    cout << "FP_FAST_FMAL: " << FP_FAST_FMAL << endl;
#endif
    cout << "FP_ILOGB0: " ; outWithHex(cout, FP_ILOGB0) << endl;
    cout << "FP_ILOGBNAN: " ; outWithHex(cout, FP_ILOGBNAN) << endl;
    cout << "HUGE_VAL: "; outWithHex(cout, HUGE_VAL) << endl;
    cout << "HUGE_VALF: "; outWithHex(cout, HUGE_VALF) << endl;
    cout << "INFINITY: "; outWithHex(cout, INFINITY) << endl;
    cout << "NAN: "; outWithHex(cout, NAN) << endl;
    cout << "quiet_NaN: "; outWithHex(cout, numeric_limits<float>::quiet_NaN()) << endl;
    cout << "signaling_NaN: "; outWithHex(cout, numeric_limits<float>::signaling_NaN()) << endl;
    cout << "denorm_min: "; outWithHex(cout, numeric_limits<float>::denorm_min()) << endl;
    cout << "lowest: "; outWithHex(cout, numeric_limits<float>::lowest()) << endl;

#ifdef __STDCPP_FLOAT16_T__
    cout << "float16_t: size=" << sizeof(float16_t) << ", MANT_DIG=" << FLT16_MANT_DIG << ", DIG=" << FLT16_DIG << ", MAX_EXP=" << FLT16_MAX_EXP << ", MIN=" ; outWithHex(cout, FLT16_MIN) << ", MAX=" ; outWithHex(cout, FLT16_MAX) << ", EPSILON=" ; outWithHex(cout, FLT16_EPSILON) << endl;
#endif
#ifdef __STDCPP_BFLOAT16_T__
    //cout << "bfloat16_t: size=" << sizeof(bfloat16_t) << ", MANT_DIG=" << BFLT16_MANT_DIG << ", DIG=" << BFLT16_DIG << ", MAX_EXP=" << BFLT16_MAX_EXP << ", MIN=" ; outWithHex(cout, BFLT16_MIN) << ", MAX=" ; outWithHex(cout, BFLT16_MAX) << ", EPSILON=" ; outWithHex(cout, BFLT16_EPSILON) << endl;
    cout << "bfloat16_t: size=" << sizeof(bfloat16_t) << ", MANT_DIG=" << numeric_limits<bfloat16_t>::digits << ", DIG=" << numeric_limits<bfloat16_t>::digits10 << ", MAX_EXP=" << numeric_limits<bfloat16_t>::max_exponent << ", MIN=" ; outWithHex(cout, numeric_limits<bfloat16_t>::min()) << ", MAX=" ; outWithHex(cout, numeric_limits<bfloat16_t>::max()) << ", EPSILON=" ; outWithHex(cout, numeric_limits<bfloat16_t>::epsilon()) << endl;
#endif
    cout << "float: size=" << sizeof(float) << ", MANT_DIG=" << FLT_MANT_DIG << ", DIG=" << FLT_DIG << ", MAX_EXP=" << FLT_MAX_EXP << ", MIN=" ; outWithHex(cout, FLT_MIN) << ", MAX=" ; outWithHex(cout, FLT_MAX) << ", EPSILON=" ; outWithHex(cout, FLT_EPSILON) << endl;
    cout << "double: size=" << sizeof(double) << ", MANT_DIG=" << DBL_MANT_DIG << ", DIG=" << DBL_DIG << ", MAX_EXP=" << DBL_MAX_EXP << ", MIN=" ; outWithHex(cout, DBL_MIN) << ", MAX=" ; outWithHex(cout, DBL_MAX) << ", EPSILON=" ; outWithHex(cout, DBL_EPSILON) << endl;
    cout << "long double: size=" << sizeof(long double) << ", MANT_DIG=" << LDBL_MANT_DIG << ", DIG=" << LDBL_DIG << ", MAX_EXP=" << LDBL_MAX_EXP << ", MIN=" ; outWithHex(cout, LDBL_MIN) << ", MAX=" ; outWithHex(cout, LDBL_MAX) << ", EPSILON=" ; outWithHex(cout, LDBL_EPSILON) << endl;
#if ALLOW_FLOAT128 && defined(__STDCPP_FLOAT128_T__)
    cout << "float128_t: size=" << sizeof(float128_t) << ", MANT_DIG=" << FLT128_MANT_DIG << ", DIG=" << FLT128_DIG << ", MAX_EXP=" << FLT128_MAX_EXP << ", MIN=" ; outWithHex(cout, FLT128_MIN) << ", MAX=" ; outWithHex(cout, FLT128_MAX) << ", EPSILON=" ; outWithHex(cout, FLT128_EPSILON) << endl;
#endif // __STDCPP_FLOAT128_T__
#if ALLOW_FLOAT128 && defined(FLT128_MANT_DIG)
    cout << "__float128: size=" << sizeof(__float128) << ", MANT_DIG=" << FLT128_MANT_DIG << ", DIG=" << FLT128_DIG << ", MAX_EXP=" << FLT128_MAX_EXP;
#ifdef FLT128_MIN
#ifdef QUADMATH_H
    // error: ambiguous overload for 'operator<<' (operand types are 'std::ostream' {aka 'std::basic_ostream<char>'} and 'const __float128')
    cout << ", MIN="; outWithHex(cout, FLT128_MIN) << ", MAX=" ; outWithHex(cout, FLT128_MAX) << ", EPSILON=" ; outWithHex(cout, FLT128_EPSILON);
#else
    cout << ", MIN="; outOnlyHex(cout, FLT128_MIN) << ", MAX=" ; outOnlyHex(cout, FLT128_MAX) << ", EPSILON=" ; outOnlyHex(cout, FLT128_EPSILON);
#endif // QUADMATH_H
#endif // FLT128_MIN
    cout << endl;
#endif // FLT128_MANT_DIG
    cout << endl;

    cout << "Run benchmark..." << endl;
    callRun();

    return 0;
}
